home *** CD-ROM | disk | FTP | other *** search
/ Power CD / Power CD ATARI-Rechner Lieben.iso / APPS / PD / EULER / PROGS / UTIL.E < prev   
Encoding:
Text File  |  1991-10-27  |  24.7 KB  |  972 lines
  1. .. The Utilities.
  2.  
  3. "Loading utilities."
  4.  
  5. Epsilon=epsilon();
  6. ResetEpsilon=epsilon();
  7. Pi=pi();
  8. E=exp(1);
  9. I=1i;
  10.  
  11. function reset
  12. ## reset() resets espilon() and some other things
  13.     global ResetEpsilon;
  14.     setepsilon(ResetEpsilon);
  15.     style(""); hold off;
  16.     color(1);
  17.     shrinkwindow();
  18.     shortformat();
  19.     return 0;
  20. endfunction
  21.  
  22. function longformat
  23. ## longformat() sets a long format for numbers
  24.     return format([20 12]);
  25. endfunction
  26.  
  27. function shortformat
  28. ## shortformat() sets a short format for numbers
  29.     return format([14 5]);
  30. endfunction;
  31.  
  32. function linspace (a,b,n)
  33. ## linspace(a,b,n) generates n+1 linear spaced points in [a,b].
  34.     if a~=b; return a; endif;
  35.     r=a:(b-a)/n:b;
  36.     r[n+1]=b;
  37.     return r;
  38. endfunction
  39.  
  40. function equispace (a,b,n)
  41. ## equispace(a,b,n) generates n+1 euqidistribution (acos) spaced values
  42. ## in [a,b].
  43.     m=(1-cos(0:pi()/n:pi()))/2;
  44.     return a+(b-a)*m;
  45. endfunction
  46.  
  47. function length (v)
  48. ## length(v) returns the length of a vector
  49.     return max(size(v));
  50. endfunction
  51.  
  52. function polydif (p)
  53. ## polydif(p) returns the polynomial p'
  54.     n=length(p);
  55.     if (n==1); return 0; endif;
  56.     return p[2:n]*(1:n-1);
  57. endfunction
  58.  
  59. function writeform (x)
  60.     if isreal(x); return printf("%25.16e",x); endif;
  61.     if iscomplex(x);
  62.         return printf("%25.16e",re(x))|printf("+%25.16ei",im(x));
  63.     endif;
  64.     if isstring(x); return x; endif;
  65.     error("Cannot print this!");
  66. endfunction
  67.  
  68. function write (x,s)
  69.     if argn()==1; s=name(x); endif;
  70.     si=size(x);
  71.     if max(si)==1; s|" = "|writeform(x)|";", return 0; endif;
  72.     s|" = [ .."
  73.     for i=1 to si[1];
  74.         for j=1 to si[2]-1;
  75.             writeform(x[i,j])|",",
  76.         end;
  77.         if i==si[1]; writeform(x[i,si[2]])|"];",
  78.         else; writeform(x[i,si[2]])|";",
  79.         endif;
  80.     end;
  81.     return 0
  82. endfunction
  83.  
  84. .. ### plot things ###
  85.  
  86. function title (text)
  87. ## title(text) plots a title to the grafik window.
  88.     ctext(text,[512 0]);
  89.     return text;
  90. endfunction
  91.  
  92. function textheight
  93. ## textheight() returns the height of a letter.
  94.     h=textsize();
  95.     return h[2];
  96. endfunction
  97.  
  98. function textwidth
  99. ## textwidth() returns the width of a letter.
  100.     h=textsize();
  101.     return h[1];
  102. endfunction
  103.  
  104. function fullwindow
  105. ## fullwindow() takes the full size (besides a title) for the
  106. ## plots.
  107.     h=textsize();
  108.     w=window();
  109.     window([12,textheight()*1.5,1011,1011]);
  110.     return w;
  111. endfunction
  112.  
  113. function shrinkwindow
  114. ## shrinkwindow() shrinks the window to allow labels.
  115.     h=textheight(); b=textwidth();
  116.     w=window();
  117.     window([6*b,1.5*h,1023-6*b,1023-3*h]);
  118.     return w;
  119. endfunction
  120.  
  121. function setplot
  122. ## setplot([xmin xmax ymin ymax]) sets the plot coordinates and holds
  123. ## the plot. Also setplot(xmin,xmax,ymin,ymax).
  124. ## setplot() resets it. 
  125.     if argn()==4; return setplot([arg1,arg2,arg3,arg4]);
  126.     else;
  127.         if argn()==0; scaling(1); 
  128.         else; error("Illegal arguments!");
  129.         endif;
  130.     endif;
  131.     return plot();
  132. endfunction
  133.  
  134. function ticks (a,b)
  135. ## returns the proper ticks to be used for intervals [a,b] and
  136. ## the factor of the ticks.
  137.     tick=10^floor(log(b-a)/log(10)-0.4);
  138.     if b-a>10*tick; tick1=tick*2; else; tick1=tick; endif;
  139.     return {(floor(a/tick1)+1)*tick1:tick1:(ceil(b/tick1)-1)*tick1,tick}
  140. endfunction
  141.  
  142. function xplot (x,y)
  143. ## xplot(x,y) or xplot(y) works like plot, but shows axis ticks.
  144.     if argn()>0;
  145.         if argn()==1;
  146.             if iscomplex(x); y=im(x); xh=re(x);
  147.             else y=x; xh=1:length(y);
  148.             endif;
  149.         else; xh=x;
  150.         endif;
  151.         p=plot(xh,y);
  152.     else
  153.         if !holding(); clg; endif;
  154.         p=plot(); frame();
  155.     endif;
  156.     {t,f}=ticks(p[1],p[2]); xgrid(t,f);
  157.     {t,f}=ticks(p[3],p[4]); ygrid(t,f);
  158.     return p;
  159. endfunction
  160.  
  161. function xmark (x,y)
  162. ## xmark(x,y) or xmark(y) works like plot, but shows axis ticks.
  163.     if argn()==1;
  164.         if iscomplex(x); y=im(x); xh=re(x);
  165.         else; y=x; xh=1:length(y);
  166.         endif;
  167.     else; xh=x;
  168.     endif;
  169.     p=mark(xh,y);
  170.     {t,f}=ticks(p[1],p[2]); xgrid(t,f);
  171.     {t,f}=ticks(p[3],p[4]); ygrid(t,f);
  172.     return p;
  173. endfunction
  174.  
  175. function fplot (ffunction,a,b)
  176. ## fplot("f",a,b,...) plots the function f in [a,b].
  177. ## The arguments ... are passed to f.
  178.     t=linspace(a,b,200);
  179.     s=ffunction(t,args(4));
  180.     return xplot(t,s)
  181. endfunction
  182.  
  183. function printscale (x)
  184.     if (abs(x)>10000) || (abs(x)<0.00001);
  185.         return printf("%12.5e",x);
  186.     else
  187.         return printf("%10.5f",x);
  188.     endif;
  189. endfunction
  190.  
  191. function xgrid
  192. ## xgrid([x0,x1,...]) draws vertical grid lines on the plot window at
  193. ## x0,x1,...
  194. ## xgrid([x0,x1,...],f) additionally writes x0/f to the axis.
  195.     c=plot(); n=length(arg1); s=scaling(0); h=holding(1);
  196.     w=window();
  197.     style("."); color(3);
  198.     ht=textheight();
  199.     loop 1 to n;
  200.         x=arg1[index()];
  201.         if (x<=c[2])&&(x>=c[1]); 
  202.             plot([x,x],[c[3],c[4]]);
  203.             if (argn()==2);
  204.                 col=w[1]+(x-c[1])/(c[2]-c[1])*(w[3]-w[1]);
  205.                 ctext(printf("%4.0lf",x/arg2),[col,w[4]+0.2*ht]);
  206.             endif;
  207.         endif;
  208.     end;
  209.     if (argn()==2);
  210.         ctext("* "|printscale(arg2),[(w[1]+w[3])/2,w[4]+1.5*ht]);
  211.     endif;
  212.     style(""); color(1); holding(h); scaling(s);
  213.     return 0;
  214. endfunction
  215.  
  216. function ygrid
  217. ## xgrid([x0,x1,...]) draws horizontal grid lines on the plot window at
  218. ## x0,x1,...
  219. ## xgrid([x0,x1,...],f) additionally writes x0/f to the axis.
  220.     c=plot(); n=length(arg1); s=scaling(0); h=holding(1);
  221.     style("."); color(3);
  222.     w=window(); wd=textwidth(); ht=textheight();
  223.     loop 1 to n;
  224.         y=arg1[index()];
  225.         if (y>=c[3])&&(y<=c[4]);
  226.             if (argn()==2);
  227.                 row=w[4]-(y-c[3])/(c[4]-c[3])*(w[4]-w[2]);
  228.                 text(printf("%4.0lf",y/arg2),[w[1]-6*wd,row-ht/2]);
  229.             endif;
  230.             plot([c[1],c[2]],[y,y]);
  231.         endif;
  232.     end;
  233.     if (argn()==2);
  234.         text("* "|printscale(arg2),[w[1]-6*wd,0]);
  235.     endif;
  236.     style(""); color(1); holding(h); scaling(s);
  237.     return 0;
  238. endfunction
  239.  
  240. function plot (x,y)
  241. ## plot(x,y) plots the values (x(i),y(i)) with the current style.
  242. ## If x is a matrix, y must be a matrix of the same size.
  243. ## The plot is then drawn for all rows of x and y.
  244. ## The plot is scaled automatically, unless hold is on.
  245. ## plot(x,y) and plot() return [x1,x2,y1,y2], where [x1,x2] is the range
  246. ## of the x-values and [y1,y2] of the y-values.
  247. ## plot(x) is the same as plot(1:length(x),x).
  248.     if argn()==1;
  249.         if iscomplex(x); return plot(re(x),im(x));
  250.         else return plot(1:length(x),x);
  251.         endif;
  252.     endif;
  253.     error("Illegal argument number!"),
  254. endfunction
  255.  
  256. function mark (x,y)
  257. ## mark(x,y) plots markers at (x(i),y(i)) according the the actual style.
  258. ## Works like plot.
  259.     if argn()==1;
  260.         if iscomplex(x); return mark(re(x),im(x));
  261.         else return mark(1:length(x),x);
  262.         endif;
  263.     endif;
  264.     error("Illegal argument number!"),
  265. endfunction
  266.  
  267. function cplot (z)
  268. ## cplot(z) plots a grid of complex numbers.
  269.     plot(re(z),im(z)); s=scaling(0); h=holding(1);
  270.     w=z'; plot(re(w),im(w)); holding(h); scaling(s);
  271.     return plot();
  272. endfunction
  273.  
  274. function plotwindow
  275. ## plotwindow() sets the plot window to the screen coordinates.
  276.     w=window();
  277.     setplot(w[1],w[3],w[2],w[4]);
  278.     return plot()
  279. endfunction
  280.  
  281. function getframe (x,y,z)
  282. ## gets the a box around all points in (x,y,z).
  283.     ex=extrema(x); ey=extrema(y); ez=extrema(z);
  284.     return [min(ex[:,1]'),max(ex[:,3]'),
  285.         min(ey[:,1]'),max(ey[:,3]'),min(ez[:,1]'),max(ez[:,3]')]
  286. endfunction
  287.  
  288. function framez0 (f)
  289.     wire([f[1],f[2],f[2],f[1],f[1]], ..
  290.         [f[3],f[3],f[4],f[4],f[3]],dup(f[5],5)');
  291.     return 0    
  292. endfunction
  293.  
  294. function framez1 (f)
  295.     wire([f[1],f[2],f[2],f[1],f[1]], ..
  296.         [f[3],f[3],f[4],f[4],f[3]],dup(f[6],5)');
  297.     return 0
  298. endfunction
  299.  
  300. function framexpyp (f)
  301.     wire([f[2],f[2]],[f[4],f[4]],[f[5],f[6]]);
  302.     return 0
  303. endfunction
  304.  
  305. function framexpym (f)
  306.     wire([f[2],f[2]],[f[3],f[3]],[f[5],f[6]]);
  307.     return 0
  308. endfunction
  309.  
  310. function framexmyp (f)
  311.     wire([f[1],f[1]],[f[4],f[4]],[f[5],f[6]]);
  312.     return 0
  313. endfunction
  314.  
  315. function framexmym (f)
  316.     wire([f[1],f[1]],[f[3],f[3]],[f[5],f[6]]);
  317.     return 0
  318. endfunction
  319.  
  320. function frame1 (f)
  321. ## draws the back part of the box (considering view)
  322.     v=view();
  323.     if v[4]>0; framez0(f);
  324.     else; framez1(f); endif;
  325.     if sin(v[3])>0;
  326.         framexmyp(f); framexmym(f);
  327.         if cos(v[3])>0; framexpyp(f); else; framexpym(f); endif;
  328.     else
  329.         framexpyp(f); framexpym(f);
  330.         if cos(v[3])>0; framexmyp(f); else; framexmym(f); endif;
  331.     endif;
  332.     return 0
  333. endfunction
  334.  
  335. function frame2 (f)
  336. ## draws the back part of the box (considering view).
  337.     v=view();
  338.     if v[4]>0; framez1(f);
  339.     else; framez0(f); endif;
  340.     if sin(v[3])>0;
  341.         if cos(v[3])>0; framexpym(f); else; framexpyp(f); endif;
  342.     else
  343.         if cos(v[3])>0; framexmym(f); else; framexmyp(f); endif;
  344.     endif;
  345.     return 0
  346. endfunction
  347.  
  348. function scaleframe (x,y,z,f,m)
  349.     s=max([f[2]-f[1],f[4]-f[3],f[6]-f[5]]);
  350.     xm=(f[2]+f[1])/2;
  351.     ym=(f[4]+f[3])/2;
  352.     zm=(f[6]+f[5])/2;
  353.     ff=m/s*2;
  354.     f=[f[1]-xm,f[2]-xm,f[3]-ym,f[4]-ym,f[5]-zm,f[6]-zm]*ff;
  355.     return {(x-xm)*ff,(y-ym)*ff,(z-zm)*ff}
  356. endfunction
  357.  
  358. function framedsolid (x,y,z,f)
  359. ## works like solid and draws a frame around the plot
  360. ## if f is specified, then the plot is scaled to fit into a cube of
  361. ## side length 2f centered at 0
  362.     frame=getframe(x,y,z);
  363.     if !holding(); clg; endif;
  364.     h=holding(1);
  365.     if argn()==4; {x1,y1,z1}=scaleframe(x,y,z,frame,f);
  366.     else; {x1,y1,z1}={x,y,z}; endif;
  367.     color(2); frame1(frame);
  368.     color(1); solid(x1,y1,z1);
  369.     color(2); frame2(frame);
  370.     color(1);
  371.     holding(h);
  372.     return frame
  373. endfunction
  374.  
  375. function framedwire (x,y,z,f)
  376. ## works like framedsolid
  377.     frame=getframe(x,y,z);
  378.     if !holding(); clg; endif;
  379.     h=holding(1);
  380.     if argn()==4; {x1,y1,z1}=scaleframe(x,y,z,frame,f);
  381.     else; {x1,y1,z1}={x,y,z}; endif;
  382.     color(2); frame1(frame);
  383.     color(1); wire(x1,y1,z1);
  384.     color(2); frame2(frame);
  385.     color(1);
  386.     holding(h);
  387.     return frame
  388. endfunction
  389.  
  390. .. ### linear algebra things ###
  391.  
  392. function id (n)
  393. ## id(n) creates a nxn identity matrix.
  394.     return diag([n n],0,1);
  395. endfunction
  396.  
  397. function inv (a)
  398. ## inv(a) produces the inverse of a matrix.
  399.     return a\id(length(a));
  400. endfunction
  401.  
  402. function fit (a,b)
  403. ## fit(a,b) computes x such that ||a.x-b||_2 is minimal.
  404.     A=conj(a');
  405.     return symmult(A,a)\(A.b')
  406. endfunction
  407.  
  408. function kernel (a)
  409. ## kernel(a) computes the kernel of the quadratic matrix a.
  410.     {aa,r,c,d}=lu(a);
  411.     n=size(a);
  412.     if length(r)==n[2] return zeros([1,n[2]])'; endif;
  413.     c1=nonzeros(c); c2=nonzeros(!c);
  414.     b=lusolve(aa[r,c1],a[r,c2]);
  415.     m=size(b);
  416.     x=zeros([n[2] m[2]]);
  417.     x[c1,:]=-b;
  418.     x[c2,:]=id(length(c2));
  419.     return x
  420. endfunction
  421.  
  422. function image (a)
  423. ## image(a) computes the image of the quadratic matrix a.
  424.     {aa,r,c,d}=lu(a);
  425.     return a[:,nonzeros(c));
  426. endfunction
  427.  
  428. function det (a)
  429. ## det(A) returns the determinant of A
  430.     {aa,r,c,d}=lu(a);
  431.     return d;
  432. endfunction
  433.  
  434. function eigenvalues (a)
  435. ## eigenvalues(A) returns the eigenvalues of A.
  436.     return polysolve(charpoly(a));
  437. endfunction
  438.  
  439. function eigenspace (a,l)
  440. ## eigenspace(A,l) returns the eigenspace of A to the eigenvaue l.
  441.     k=kernel(a-l*id(length(a)));
  442.     if k==0; error("No eigenvalue found!"); endif;
  443.     si=size(k);
  444.     loop 1 to si[2];
  445.         x=k[:,index()];
  446.         k[:,index()]=x/sqrt(x'.x);
  447.     end;
  448.     return k;
  449. endfunction
  450.  
  451. function eigen (A)
  452. ## eigen(A) returns the eigenvectors and a basis of eigenvectors of A.
  453. ## {l,x}=eigen(A), where l is a vector of distinct eigenvalues and
  454. ## x is a basis of eigenvectors.
  455.     l=eigenvalues(A);
  456.     s=eigenspace(A,l[1]);
  457.     si=size(s); v=l[1]; l=eigenremove(l,si[2]);
  458.     repeat;
  459.         if min(size(l))==0; break; endif;
  460.         ll=l[1]; sp=eigenspace(A,ll);
  461.         si=size(sp); l=eigenremove(l,si[2]);
  462.         s=s|sp; v=v|ll;
  463.     end;
  464.     return {v,s}
  465. endfunction
  466.  
  467.  
  468. function eigenspace1
  469. ## eigenspace1(A,l) returns the eigenspace of A to the eigenvalue l.
  470. ## returns {x,l1}, where l1 should be an improvement over l, and
  471. ## x contains the eigenvectors as columns.
  472.     eps=epsilon();
  473.     repeat;
  474.         k=kernel(arg1-arg2*id(length(arg1)));
  475.         if k==0; else; break; endif;
  476.         if epsilon()>1 break; endif;
  477.         setepsilon(100*epsilon());
  478.     end;
  479.     if k==0; error("No eigenvalue found!"); endif;
  480.     setepsilon(eps);
  481.     {dummy,l}=eigennewton(arg1,k[:,1],arg2);
  482.     eps=epsilon();
  483.     repeat;
  484.         k=kernel(arg1-l*id(length(arg1)));
  485.         if k==0; else; break; endif;
  486.         if epsilon()>1 break; endif;
  487.         setepsilon(100*epsilon());
  488.     end;
  489.     if k==0; error("No eigenvalue found!"); endif;
  490.     setepsilon(eps);
  491.     si=size(k);
  492.     loop 1 to si[2];
  493.         x=k[:,index()];
  494.         k[:,index()]=x/sqrt(x'.x);
  495.     end;
  496.     return {k,l};
  497. endfunction
  498.  
  499. function eigenremove
  500. ## helping function.
  501.     {l,n}={arg1,arg2};
  502.     k=length(l); l1=l[2:k];
  503.     if n==1; return l1; endif;
  504.     v=abs(l[1]-l1);
  505.     {vs,i}=sort(v);
  506.     return(l1[i[n:k]]);
  507. endfunction
  508.  
  509. function eigen1
  510. ## eigen1(A) returns the eigenvectors and a basis of eigenvectors of A.
  511. ## {l,x}=eigen(A), where l is a vector of distinct eigenvalues and
  512. ## x is a basis of eigenvectors. Improved version of eigen.
  513.     l=eigenvalues(arg1);
  514.     {s,ll}=eigenspace1(arg1,l[1]); l[1]=ll;
  515.     si=size(s); v=l[1]; l=eigenremove(l,si[2]);
  516.     repeat;
  517.         if min(size(l))==0; break; endif;
  518.         ll=l[1];
  519.         {sp,ll}=eigenspace1(arg1,ll); l[1]=ll;
  520.         si=size(sp); l=eigenremove(l,si[2]);
  521.         s=s|sp; v=v|ll;
  522.     end;
  523.     return {v,s}
  524. endfunction
  525.  
  526. function eigennewton
  527. ## eigennewton(a,x,l) does a newton step to improve the eigenvalue l
  528. ## of a and the eigenvector x.
  529. ## returns {x1,l1}.
  530.     a=arg1; x=arg2; x=x/sqrt(x'.x); n=length(a);
  531.     d=((a-arg3*id(n))|-x)_(2*x'|0);
  532.     b=d\((a.x-arg3*x)_0);
  533.     return {x-b[1:n],arg3-b[n+1]};
  534. endfunction
  535.  
  536. hilbfactor=5*7*8*27*11*13*17*19*23*29;
  537.  
  538. function hilb (n)
  539. ## hilb(n) produces the nxn-Hilbert matrix, multiplied by hilbfactor.
  540.     global hilbfactor;
  541.     {i,j}=field(1:n,1:n);
  542.     return hilbfactor/(i+j-1);
  543. endfunction
  544.  
  545. .. ### polynomial fit ##
  546.  
  547. function polyfit (xx,yy,n)
  548. ## fit(x,y,degree) fits a polynomial in L_2-norm to (x,y).
  549.     A=ones(size(xx))_dup(xx,n); A=cumprod(A');
  550.     return fit(A,yy)';
  551. endfunction
  552.  
  553. function field (x,y)
  554. ## field(x,y) returns {X,Y} such that the X+i*Y is a rectangular
  555. ## grid in C containing the points x(n)+i*y(m).
  556.     return {dup(x,length(y)),dup(y',length(x))};
  557. endfunction
  558.  
  559. function totalsum (A)
  560. ## totalsum(a) computes the sum of all elements of a
  561.     return sum(sum(A)');
  562. endfunction
  563.  
  564. function sinh
  565. ## sinh(x) = (exp(x)-exp(-x))/2
  566.     h=exp(arg1);
  567.     return (h-1/h)/2;
  568. endfunction
  569.  
  570. function cosh
  571. ## cosh(x) = (exp(x)+exp(-x))/2
  572.     h=exp(arg1);
  573.     return (h+1/h)/2;
  574. endfunction
  575.  
  576. function arsinh
  577. ## arsinh(z) = log(z+sqrt(z^2+1))
  578.     return log(arg1+sqrt(arg1*arg1+1))
  579. endfunction
  580.  
  581. function arcosh
  582. ## arcosh(z) = log(z+(z^2-1))
  583.     return log(arg1+sqrt(arg1*arg1-1))
  584. endfunction
  585.  
  586. function heun (ffunction,t,y0)
  587. ## y=heun("f",t,y0,...) solves the differential equation y'=f(t,y).
  588. ## f(t,y,...) must be a function. 
  589. ## y0 is the starting value.
  590.     n=length(t);
  591.     y=dup(y0,n);
  592.     loop 1 to n-1;
  593.         h=t[#+1]-t[#];
  594.         yh=y[#]; xh=t[#];
  595.         k1=ffunction(xh,yh,args(4));
  596.         k2=ffunction(xh+h/2,yh+h/2*k1,args(4));
  597.         k3=ffunction(xh+h,yh+2*h*k2-h*k1,args(4));
  598.         y[#+1]=yh+h/6*(k1+4*k2+k3);
  599.     end;
  600.     return y';
  601. endfunction
  602.  
  603. solveeps=epsilon();
  604.  
  605. function bisect (ffunction,a,b)
  606. ## bisect("f",a,b,...) uses the bisection method to find a root of 
  607. ## f in [a,b].
  608.     global solveeps;
  609.     if ffunction(b,args(4))<0;
  610.         b1=a; a1=b;
  611.     else;
  612.         a1=a; b1=b;
  613.     endif;
  614.     if ffunction(b1,args(4))<0 error("No zero in interval"); endif;
  615.     repeat
  616.         m=(a1+b1)/2;
  617.         if abs(a1-b1)<solveeps; return m; endif;
  618.         if ffunction(m,args(4))>0; b1=m; else a1=m; endif;
  619.     end;
  620. endfunction
  621.  
  622. function secant (ffunction,a,b)
  623. ## secant("f",a,b,...) uses the secant method to find a root of f in [a,b]
  624.     global solveeps;
  625.     x0=a; x1=b; y0=ffunction(x0,args(4)); y1=ffunction(x1,args(4));
  626.     repeat
  627.         x2=x1-y1*(x1-x0)/(y1-y0);
  628.         if abs(x2-x1)<solveeps; break; endif;
  629.         x0=x1; y0=y1; x1=x2; y1=ffunction(x2,args(4));
  630.     end;
  631.     return x2
  632. endfunction
  633.  
  634. function simpson (ffunction,a,b)
  635. ## simpson("f",a,b) or simpson("f",a,b,n,...) integrates f in [a,b] with
  636. ## the simpson method. f must be able to evaluate a vector.
  637.     if argn()>=4; n=arg4; else; n=50; endif;
  638.     t=linspace(a,b,2*n);
  639.     s=ffunction(t,args(5));
  640.     ff=4-mod(1:2*n+1,2)*2; ff[1]=1; ff[2*n+1]=1;
  641.     return sum(ff*s)/3*(t[2]-t[1]);
  642. endfunction
  643.  
  644. rombergeps=epsilon();
  645.  
  646. function romberg(ffunction,a,b,m)
  647. ## romberg(f,a,b) computes the Romberg integral of f in [a,b].
  648. ## romberg(f,a,b,m,...) specifies h=(b-a)/m/2^k for k=1,...
  649. ## ... is passed to f(x,...)
  650.     global rombergeps;
  651.     if argn()==3; m=10; endif;
  652.     y=ffunction(linspace(a,b,m),args(5)); h=(b-a)/m;
  653.     y[1]=y[1]/2; y[m+1]=y[m+1]/2; I=sum(y);
  654.     S=I*h; H=h^2; Intalt=S;
  655.     repeat;
  656.         I=I+sum(ffunction(a+h/2:h:b,args(5))); h=h/2;
  657.         S=S|I*h;
  658.         H=H|h^2;
  659.         Int=interpval(H,interp(H,S),0);
  660.         if abs(Int-Intalt)<rombergeps; break; endif;
  661.         Intalt=Int;
  662.     end;
  663.     return Intalt
  664. endfunction
  665.  
  666. gaussz = [ ..
  667. -9.7390652851717172e-0001,
  668. -8.6506336668898451e-0001,
  669. -6.7940956829902441e-0001,
  670. -4.3339539412924719e-0001,
  671. -1.4887433898163121e-0001,
  672.  1.4887433898163121e-0001,
  673.  4.3339539412924719e-0001,
  674.  6.7940956829902440e-0001,
  675.  8.6506336668898451e-0001,
  676.  9.7390652851717172e-0001];
  677. gaussa = [ ..
  678.  6.6671344308688139e-0002,
  679.  1.4945134915058059e-0001,
  680.  2.1908636251598205e-0001,
  681.  2.6926671930999635e-0001,
  682.  2.9552422471475288e-0001,
  683.  2.9552422471475287e-0001,
  684.  2.6926671930999635e-0001,
  685.  2.1908636251598205e-0001,
  686.  1.4945134915058059e-0001,
  687.  6.6671344308688137e-0002];
  688.  
  689. function gauss10 (ffunction,a,b)
  690. ## gauss10("f",a,b,...)
  691. ## evaluates the gauss integration with 10 knots an [a,b]
  692.     global gaussa,gaussz;
  693.     z=a+(gaussz+1)*(b-a)/2;
  694.     return sum(gaussa*ffunction(z,args(4)))*(b-a)/2;
  695. endfunction
  696.  
  697. function gauss (ffunction,a,b,n)
  698. ## gauss("f",a,b) gauss integration with 10 knots
  699. ## gauss("f",a,b,n,...) subdivision into n subintervals.
  700. ## a and b may be vectors.
  701.     si=size(a,b); R=zeros(si);
  702.     if argn()==3;
  703.         loop 1 to prod(si);
  704.             sum=0;
  705.             sum=sum+gauss10(ffunction,a{#},b{#});
  706.             R{#}=sum;
  707.         end;
  708.     else;
  709.         loop 1 to prod(si);
  710.             h=linspace(a{#},b{#},n); sum=0;
  711.             loop 1 to n;
  712.                 sum=sum+gauss10(ffunction,h{#},h{#+1},args(5));
  713.             end;
  714.             R{#}=sum;
  715.         end;
  716.     endif;
  717.     return R
  718. endfunction
  719.  
  720. .. ### Broyden ###
  721.  
  722. broydenmax=50;
  723.  
  724. function broyden (ffunction,xstart,A)
  725. ## broyden("f",x) or broyden("f",x,a,...) finds a zero of f.
  726. ## x is the starting value and a is a approximation of the jacobian.
  727. ## if a is 0, it is neglected.
  728. ## ... is passed to f(x,...)
  729.     global solveeps;
  730.     x=xstart; n=length(x);
  731.     global broydenmax;
  732.     delta=sqrt(solveeps);
  733.     if argn()==2; A=0; endif;
  734.     if A==0;
  735.         A=zeros([n n]);
  736.         y=ffunction(x,args(4));
  737.         loop 1 to n;
  738.             x0=x; x0[#]=x0[#]+delta;
  739.             A[:,#]=(ffunction(x0,args(4))-y)'/delta;
  740.         end;
  741.     endif;
  742.     itercount=0;
  743.     y=ffunction(x,args(4));
  744.     repeat
  745.         d=-A\y'; x=x+d';
  746.         y1=ffunction(x,args(4)); q=y1-y;
  747.         A=A+((q'-A.d).d')/(d'.d);
  748.         if (max(abs(d))<solveeps || itercount>broydenmax) break; endif;
  749.         y=y1;
  750.         itercount=itercount+1;
  751.     end;
  752.     if (itercount>broydenmax) error("Iteration failed"); endif;
  753.     return x;
  754. endfunction
  755.  
  756. function map (ffunction,t)
  757. ## map("f",t,...) evaluates the function f at all points of the vector t.
  758. ## usually this is the same as f(t,...), but sometimes functions do not
  759. ## accept vectors as input.
  760.     si=size(t); s=zeros(si);
  761.     loop 1 to prod(si); s{#}=ffunction(t{#},args(3)); end;
  762.     return s;
  763. endfunction
  764.  
  765. itereps=epsilon();
  766.  
  767. function iterate (ffunction,x,n)
  768. ## iterate("f",x,n,...) iterate the function f(x,...) n times,
  769. ## starting with the point x.
  770. ## if n is missing or n is 0, then the iteration stops at a fixed point.
  771. ## Returns the value after n iterations, and its history.
  772.     global itereps;
  773.     if argn()<3; n=0; endif;
  774.     if n==0;
  775.         R=x; Rold=x;
  776.         repeat
  777.             Rnew=ffunction(Rold,args(4));
  778.             R=R_Rnew;
  779.             if max(abs(Rold-Rnew))<itereps; return Rnew; endif;
  780.             Rold=Rnew;
  781.         end;
  782.     else;
  783.         R=zeros(n,length(x));
  784.         R[1]=x;
  785.         loop 2 to n; R[#]=ffunction(R[#-1],args(4)); end;
  786.         return {R[n],R'}
  787.     endif;
  788. endfunction    
  789.  
  790. .. ### Remez algorithm ###
  791.  
  792. function remezhelp (x,y,n)
  793. ## {t,d,h}=remezhelp(x,y,deg) is the case, where x are exactly deg+2
  794. ## points.
  795.     d1=interp(x,y); d2=interp(x,mod(1:n+2,2)*2-1);
  796.     h=d1[n+2]/d2[n+2];
  797.     d=d1-d2*h;
  798.     return {x[1:n+1],d[1:n+1],h};
  799. endfunction
  800.  
  801. remezeps=0.00001;
  802.  
  803. function remez
  804. ## {t,d,h,r}=remez(x,y,deg) computes the divided difference form
  805. ## of the polynomial best approximation to y at x of degree deg.
  806. ## remez(x,y,deg,1) does the same thing with tracing.
  807. ## h is the discrete error (with sign), and x the rightmost point,
  808. ## which is missing in t.
  809.     global remezeps;
  810.     if argn()==4; trace=1; else; trace=0; endif;
  811.     x=arg1; y=arg2; deg=arg3; n=length(x);
  812.     ind=linspace(1,n,deg+1); h=0;
  813.     repeat
  814.         {t,d,hnew}=remezhelp(x[ind],y[ind],deg);
  815.         r=y-interpval(t,d,x); hh=hnew;
  816.         if trace;
  817.             xind=x[ind];
  818.             plot(x,r); xgrid(xind); ygrid([-h,0,h]);
  819.             title("Weiter mit Taste"); wait(180);
  820.         endif;
  821.         indnew=ind; rr=ind;
  822.         rr[1]=r[ind[1]];
  823.         rr[2]=r[ind[deg+2]];
  824.         loop 2 to deg+1;
  825.             e=extrema(r[ind[index()-1]:ind[index()+1]]);
  826.             if (hh>0);
  827.                 indnew[index()]=e[2]+ind[index()-1]; rr[index()]=e[1];
  828.             else
  829.                 indnew[index()]=e[4]+ind[index()-1]; rr[index()]=e[3];
  830.             endif;
  831.             hh=-hh;
  832.         end;
  833.         ind=indnew;
  834.         if (abs(hnew)<(1+remezeps)*abs(h)) break; endif;
  835.         h=hnew;
  836.     end;
  837.     {t,d,h}=remezhelp(x[ind],y[ind],deg);
  838.     if trace;
  839.         xind=x[ind];
  840.         plot(x,y-interpval(t,d,x)); xgrid(xind); ygrid([-h,0,h]);
  841.         title("Weiter mit Taste"); wait(180);
  842.     endif;
  843.     return {t,d,h,x[ind[deg+2]]};
  844. endfunction
  845.  
  846. .. ### Natural spline ###
  847.  
  848. function spline (x,y)
  849. ## spline(x,y) defines the natural Spline at points x(i) with
  850. ## values y(i). It returns the second derivatives at these points.
  851.     n=length(x);
  852.     h=x[2:n]-x[1:n-1];
  853.     s=h[2:n-1]+h[1:n-2];
  854.     l=h[2:n-1]/s;
  855.     A=diag([n,n],0,2);
  856.     A=setdiag(A,1,0|l);
  857.     A=setdiag(A,-1,1-l|0);
  858.     b=6/s*((y[3:n]-y[2:n-1])/h[2:n-1]-(y[2:n-1]-y[1:n-2])/h[1:n-2]);
  859.     return (A\(0|b|0)')';
  860. endfunction
  861.  
  862. function splineval (x,y,h,t)
  863. ## splineval(x,y,h,t) evaluates the interpolation spline for
  864. ## (x(i),y(i)) with second derivatives h(i) at t.
  865.     p1=find(x,t); p2=p1+1;
  866.     y1=y[p1]; y2=y[p2];
  867.     x1=x[p1]; x2=x[p2]; f=(x2-x1)^2;
  868.     h1=h[p1]*f; h2=h[p2]*f;
  869.     b=h1/2; c=(h2-h1)/6;
  870.     a=y2-y1-b-c;
  871.     d=(t-x1)/(x2-x1);
  872.     return y1+d*(a+d*(b+c*d));
  873. endfunction
  874.  
  875. .. ### use zeros,... the usual way ###
  876.  
  877. function ctext
  878. ## ctext(s,[c,r]) plots the centered string s at the coordinates (c,r).
  879. ## Also ctext(s,c,r).
  880.     if argn()==3; return ctext(arg1,[arg2,arg3]); endif;
  881.     error("Illegal argument number!"),
  882. endfunction
  883.  
  884. function text
  885. ## text(s,[c,r]) plots the centered string s at the coordinates (c,r).
  886. ## Also ctext(s,c,r).
  887.     if argn()==3; return text(arg1,[arg2,arg3]); endif;
  888.     error("Illegal argument number!"),
  889. endfunction
  890.  
  891. function diag
  892. ## diag([n,m],k,v) returns a nxm matrix A, containing v on its k-th
  893. ## diagonal. v may be a vector or a real number. Also diag(n,m,k,v);
  894. ## diag(A,k) returns the k-th diagonal of A.
  895.     if argn()==4; return diag([arg1,arg2],arg3,arg4); endif;
  896.     error("Illegal argument number!"),
  897. endfunction
  898.  
  899. function format
  900. ## format([n,m]) sets the number output format to m digits and a total
  901. ## width of n. Also format(n,m);
  902.     if argn()==2; return format([arg1,arg2]); endif;
  903.     error("Illegal argument number!"),
  904. endfunction
  905.  
  906. function normal
  907. ## normal([n,m]) returns a nxm matrix of unit normal distributed random
  908. ## values. Also normal(n,m).
  909.     if argn()==2; return normal([arg1,arg2]); endif;
  910.     error("Illegal argument number!"),
  911. endfunction
  912.  
  913. function random
  914. ## random([n,m]) returns a nxm matrix of uniformly distributed random 
  915. ## values in [0,1]. Also random(n,m).
  916.     if argn()==2; return random([arg1,arg2]); endif;
  917.     error("Illegal argument number!"),
  918. endfunction
  919.  
  920. function ones
  921. ## ones([n,m]) returns a nxm matrix with all elements set to 1.
  922. ## Also ones(n,m).
  923.     if argn()==2; return ones([arg1,arg2]); endif;
  924.     error("Illegal argument number!"),
  925. endfunction
  926.  
  927. function zeros
  928. ## zeros([n,m]) returns a nxm matrix with all elements set to 0.
  929. ## Also zeros(n,m).
  930.     if argn()==2; return zeros([arg1,arg2]); endif;
  931.     error("Illegal argument number!"),
  932. endfunction
  933.  
  934. function matrix
  935. ## matrix([n,m],x) returns a nxm matrix with all elements set to x.
  936. ## Also matrix(n,m,x).
  937.     if argn()==3; return matrix([arg1,arg2],arg3); endif;
  938.     error("Illegal argument number!"),
  939. endfunction
  940.  
  941. function view
  942. ## view([distance, tele, angle1, angle2]) sets the perspective for
  943. ## solid and view. distance is the eye distance, tele a zooming factor.
  944. ## angle1 is the angle from the negativ y-axis to the positive x-axis.
  945. ## angle2 is the angle to the positive z-axis (the height of the eye).
  946. ## Also view(d,t,a1,a2).
  947. ## view() returns the values of view.
  948.     if argn()==4; return view([arg1,arg2,arg3,arg4]); endif;
  949.     error("Illegal argument number!"),
  950. endfunction
  951.  
  952. function window
  953. ## window([c1,r1,c2,r2]) sets a plotting window. The cooridnates must
  954. ## be screen coordimates. Also window(c1,r1,c2,r2).
  955. ## window() returns the active window coordinates.
  956.     if argn()==4; return window([arg1,arg2,arg3,arg4]); endif;
  957.     error("Illegal argument number!"),
  958. endfunction
  959.  
  960. .. *** directory ***
  961.  
  962. function dir(pattern)
  963. ## dir(pattern) displays a directory.
  964. ## dir() is the same as dir("*.*").
  965.     if argn()==0; pattern="*.*"; endif;
  966.     s=searchfile(pattern), if stringcompare(s,"")==0; return 0; endif;
  967.     repeat
  968.         s=searchfile(), if stringcompare(s,"")==0; return 0; endif;
  969.     end;
  970. endfunction
  971.  
  972.